{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noise filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sdm as sdmlib\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image, ImageDraw, ImageFont\n",
    "import urllib, cStringIO\n",
    "import random\n",
    "from IPython.core.display import display, Image as IPythonImage\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "width = 30\n",
    "height = 30\n",
    "noise_flip = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_img(letter='A'):\n",
    "    img = Image.new('RGBA', (30, 30), (255, 255, 255))\n",
    "    font = ImageFont.truetype('Arial.ttf', 30)\n",
    "    draw = ImageDraw.Draw(img)\n",
    "    draw.text((5, 0), letter, (0, 0, 0), font=font)\n",
    "    return img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_noise_add(img, p=0.15, flip=False):\n",
    "    img2 = img.copy()\n",
    "    draw = ImageDraw.Draw(img2)\n",
    "    for py in xrange(height):\n",
    "        for px in xrange(width):\n",
    "            if random.random() < p:\n",
    "                if flip:\n",
    "                    pixel = img.getpixel((px, py))\n",
    "                    value = sum([int(x/255+0.5) for x in pixel[:3]])//3\n",
    "                    assert value == 0 or value == 1\n",
    "                    value = (1 - value)*255\n",
    "                    draw.point((px, py), fill=(value, value, value))\n",
    "                else:\n",
    "                    draw.point((px, py), fill=(0, 0, 0))\n",
    "    return img2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAC7CAYAAAB1qmWGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEZtJREFUeJzt3VuIXVWex/Hfv+2Yh9GHaKqT4GUS2hANSCfU0RnoRmMyDl7A6It0SBoHAmlkFMV4iYLpRhxR7HbmoYeGKgypiGPTUl7yEGaUXDANg3gqEW9REyRGi1wq+GANgcSk/vNQx6bMWtW16+x99tl71fcDoar+tc9Za1f982fX+Z+1trm7AAD196NuTwAAUAwKOgAkgoIOAImgoANAIijoAJAICjoAJIKCDgCJoKADQCJyFXQzu8XMPjOzQ2a2qahJAd1GbqOOrN2VomZ2gaTPJd0s6WtJ70la4+6fFDc9oHzkNurqxzkee72kQ+7+hSSZ2Z8krZY0adLPnTvXFy5cmGNIYHKHDx/WyZMnrYCnqmRuDw0NReO9vb0dHTev2LzzzrkTz1nlcbPmdp6CfpmkryZ8/bWkf/hbD1i4cKGazWaOIYHJNRqNop6qkrltFv//XPX/U7F5551zJ56zyuNmze2ON0XNbIOZNc2sOTIy0unhgNKQ26iaPAV9WNIVE76+vBX7AXfvc/eGuzd6enpyDAeUhtxGLeUp6O9JWmxmi8zsQkm/lLS9mGkBXVVIbptZ8G8mcvfgXxWfs4xxO50Tbb+G7u5nzew+Sf8j6QJJW9z948JmBnQJuY26ytMUlbvvkLSjoLkAlUFuo45YKQoAiaCgA0Aicr3kAmByRTfqOtH4izXluM9w58R+tkX+DrhCB4BEUNABIBEUdABIBAUdABJBUxQowNDQUNDcKqu5mKepVocGaOqN2yLPhSt0AEgEBR0AEkFBB4BEUNABIBEUdKAAvb29wbaqRW+VGns+M+vaVrJlqfr5TfZ76QYKOgAkgoIOAImgoANAIijoAJCIXCtFzeywpFFJ5ySddfdGEZMCuq2I3M7TvMu7OrJbqyuzNgOr1tjMo0rnUsTS/5vc/WQBzwNUDbmNWuElFwBIRN6C7pLeMrMhM9tQxISAiiC3UTt5X3L5hbsPm9lPJL1tZp+6+zsTD2j9Z9ggSVdeeWXO4YDSkNuonVxX6O4+3Pp4QtLrkq6PHNPn7g13b/T09OQZDigNuY06avsK3cz+TtKP3H209fk/S3qqsJkl4tSpU0Fs/vz5QWx0dDTzcz733HNB7NFHH53exDCpdnI7th96TFn7lHfrnRcXX3xxEIvl9mQ/q6y5nfoe6e1uHZDnJZd5kl5vDfxjSf/l7v+d4/mAqiC3UUttF3R3/0LSzwqcC1AJ5DbqirctAkAiKOgAkAhuEt1hg4ODQWw6DdCY/v7+IPbII48EsW7tyTwT9fb2qtlsdnsaf9WtpmHW3J5sLosXLw5isdzu5rYKeWQd+/xYo5Ft5wmu0AEgERR0AEgEBR0AEkFBB4BE0BTtsK1btxb+nIcOHQpiu3btCmKrVq0qfGyUYzqNuzKa37ExVq5cGcRic4w9djpzznpsuw3HMnV6bK7QASARFHQASAQFHQASQUEHgETQFC3QkSNHgtiePXsyPfaOO+4IYtu3b888dl9fXxCjKTozlLFq8ssvvwxiixYtyvR8eXM7q241Oydr2nZjPlyhA0AiKOgAkAgKOgAkYsqCbmZbzOyEmX00IXaJmb1tZgdbH+d0dppA8chtpCZLU3SrpD9I2jYhtknSTnd/1sw2tb5+rPjp1cvAwEAQGxsby/TY+++/P4h99dVX0WP3798fxN54440gNjIyEsS4mfEPbFVFc7ushlrWcbqZ27Nnzw5ip0+fzjR20crakrrdcaa8Qnf3dyR9c154taTvf8MDku5sa3Sgi8htpKbd19DnufvR1ufHNH5TXSAF5DZqK3dT1Mf/Zpv07zYz22BmTTNrxl4CAKqK3EbdtFvQj5vZAklqfTwx2YHu3ufuDXdv8PotaoDcRm21u1J0u6R7JD3b+vhmYTOqsW3btk19kKQ5c8I3Ttx4441BbN26ddHHxxpHZ86cCWKxrXtj92fED3Q0t8u6n2XWcbI236666qpMx+XN7X379mUa5/nnnw9iRed2N+892rF7iprZK5L+V9ISM/vazNZrPNlvNrODkv6p9TVQK+Q2UjPlFbq7r5nkW2wUglojt5EaVooCQCIo6ACQCLbPbdPevXuDWOxenzF33XVXEJs1a1YQW7t2bfTxjz0WLlw8e/ZsEOvv7w9iDz/8cBAra/Ub8jUmp9OQy3ps1uOyzjFvbseO/e6774JYt3K7m43SLLhCB4BEUNABIBEUdABIBAUdABJBU7RNsVWYWd17772Zjps3L74v1OrVq4PY4OBgEDt48GAQ2717dxBbuXJlpvmgM/I2JrM2WrOOs379+kzHxXQit7M2NmO5HbuvbtEN4yrhCh0AEkFBB4BEUNABIBEUdABIBE3RDE6dOhXEXn311UyPXbFiRRDLuhXmZDZt2hTEYk3RmL6+viBGU7QeJmvSZW0aZm2UZs3tTjQN8+R2LI/r2NjMgyt0AEgEBR0AEkFBB4BEZLlj0RYzO2FmH02I/dbMhs3s/da/2zo7TaB45DZSk6UpulXSHySdf8PMf3f33xU+owqKNWVGR0czPbYT9/CMNVVjDaFdu3YFsddffz2InTx5MojNnTu3zdnVylZ1MLc7sS1u0Y9/6aWXgljW3N6xY0cQu/3224PYdOaXJ7dnz54dxIaHh4NYyrk95RW6u78j6ZsS5gKUitxGavK8hn6fmX3Q+rM1vNU3UF/kNmqp3YL+R0k/lbRM0lFJv5/sQDPbYGZNM2uOjIy0ORxQGnIbtdVWQXf34+5+zt3HJPVLuv5vHNvn7g13b/T09LQ7T6AU5DbqrK2CbmYLJnx5l6SPJjsWqBNyG3U25btczOwVSSskzTWzryX9RtIKM1smySUdlvTrDs6x67Lufb506dIgduuttxY8m7jYkunYOwHOnDkTxGLnF7vhbmo6ndt53n1S1s2Ii87tTsxx586dmY6L/cxi5xd751kn5t2NG0pPWdDdfU0k/GIH5gKUitxGalgpCgCJoKADQCIo6ACQCPZDP8+RI0eC2J49ezI9duPGjUEs6/LvvG6++eYgFltG3Ww2g1h/f38QmwlN0W4oazuAmCrl9nRueB1z3XXXZXrskiVLMh3XiQZmN/Zi5wodABJBQQeARFDQASARFHQASARN0fMMDAwEsbGxsUyPXb9+faZY1Xz++edBbPfu3UHspptuKmM6SevmTYvrkNt5msZZHxvL7VRuJs0VOgAkgoIOAImgoANAIijoAJAImqLn2bbt/PsFz0x9fX1BjKZovT355JNBbPPmzV2YSXmKbnZ2Y0vc6eAKHQASQUEHgERQ0AEgEVMWdDO7wsx2m9knZvaxmT3Qil9iZm+b2cHWxzmdny5QHHIbqcnSFD0raaO77zOziyUNmdnbkv5F0k53f9bMNknaJOmxzk21eHv37g1ihw4d6sJMque1114LYidPngxic+fOLWM6nTKjcvuGG27owkymJ+sK0KyrQmPHXXjhhUFseHg4iMVyezorVIu+p2wWU16hu/tRd9/X+nxU0gFJl0laLen7tcQDku5sawZAl5DbSM20XkM3s4WSlkt6V9I8dz/a+tYxSfMKnRlQInIbKchc0M3sIkmDkh50928nfs/H/7aI/n1hZhvMrGlmzZGRkVyTBTqB3EYqMhV0M5ul8YR/2d2/f3H1uJktaH1/gaQTsce6e5+7N9y90dPTU8ScgcKQ20jJlE1RG391/kVJB9z9hQnf2i7pHknPtj6+2ZEZdtDWrVvbfuyBAweC2NVXX51jNuWIrRZ8+umng9iZM2eCWGz71di9JuuiKrndidWHeXK7SisfpXJWe+bJ7cnml+f3ev5xsfsDx2R5l8vPJf1K0odm9n4r9oTGk/3PZrZe0peS7s40IlAd5DaSMmVBd/e/SJrsPTSrip0OUB5yG6lhpSgAJIKCDgCJsDIbII1Gw5vNZmnjTXTq1KkgNn/+/CA2OjoaxK699tog9sEHHxQzsZJ99tlnQSxrM3fJkiVB7NNPP809p6I0Gg01m832ltjlHzvI7bK2Wi0jt9tduSjlP+esP8eUcnuSc57yl8AVOgAkgoIOAImgoANAIijoAJCIGXNP0cHBwSAWaxLF3H13OutKYs2f3t7eIDY0NBTEYk2nPXv2BLEVK1a0NbeZoBON0qy5HRsntko4Juu2sdPZXjbrOFnlye1YA7Sbud3uSlGu0AEgERR0AEgEBR0AEkFBB4BEzJimaJ7tRFNqisasW7cuiMUaRzF9fX1BjKbouLJWYWfN7VhzMtbozipPozTvOFllze3YvNesWRPEqp7bXKEDQCIo6ACQCAo6ACRiyoJuZleY2W4z+8TMPjazB1rx35rZsJm93/p3W+enCxSH3EZqptw+t3WT3AXuvs/MLpY0JOlOjd+W6//c/XdZBytj+9wjR45E44sWLQpiY2NjQWzZsmVBbP/+/fknVmHHjh0LYpdffnkQO3fuXBCbPXt2EBseHg5il156aZuzy2662+dWJbezNhKnk9ux39Xy5cuDWCy3y9r2twwzLbez3ILuqKSjrc9HzeyApMvyTxHoLnIbqZnWa+hmtlDScknvtkL3mdkHZrbFzOYUPDegNOQ2UpC5oJvZRZIGJT3o7t9K+qOkn0papvGrnN9P8rgNZtY0s+bIyEgBUwaKRW4jFZkKupnN0njCv+zur0mSux9393PuPiapX9L1sce6e5+7N9y90dPTU9S8gUKQ20hJlne5mKQXJR1w9xcmxBdMOOwuSR8VPz2gc8htpCbL0v+fS/qVpA/N7P1W7AlJa8xsmSSXdFjSrzsyw2kaGBiIxmPvaIlJfZl/TOyGwqtWrQpib731VhA7ffp0EIv9Dh566KE2Z9dRlcjtrEvon3rqqejjY7kde/wzzzyT6bg8e5qXtR1ATNZxis7tjRs3RueT5xzbvSl3lne5/EVS7Nl3tDUiUBHkNlLDSlEASAQFHQASQUEHgEQktx/6tm3bcj1+JjZFY9auXRvEYo2jmP7+/iBW0aZoZU2nobZ58+ZMx8Vy+/HHHw9ieRqgeY7LK+s4WXM79nzXXHNN2+NOBzeJBoAZjoIOAImgoANAIijoAJCIKfdDL1IZ+6Fj5prufuhFMrPgP1I3V00WrWqrPas+RtFjZ81trtABIBEUdABIBAUdABJBQQeARCS3UhToht7eXmVp+Fe9+TaZoreCnez5yvj5dLMJnWdr4iy4QgeARFDQASARFHQASAQFHQASUepKUTMbkfSlpLmSTpY2cGdxLtXx9+7e042Bye3Kq/u5ZMrtUgv6Xwc1a7p7tg1+K45zwUQp/Qw5l/rhJRcASAQFHQAS0a2C3telcTuBc8FEKf0MOZea6cpr6ACA4vGSCwAkovSCbma3mNlnZnbIzDaVPX4eZrbFzE6Y2UcTYpeY2dtmdrD1cU4355iFmV1hZrvN7BMz+9jMHmjFa3cuVUJud99Mz+1SC7qZXSDpPyXdKmmppDVmtrTMOeS0VdIt58U2Sdrp7osl7Wx9XXVnJW1096WS/lHSv7Z+D3U8l0ogtytjRud22Vfo10s65O5fuPsZSX+StLrkObTN3d+R9M154dWSBlqfD0i6s9RJtcHdj7r7vtbno5IOSLpMNTyXCiG3K2Cm53bZBf0ySV9N+PrrVqzO5rn70dbnxyTN6+ZkpsvMFkpaLuld1fxcuozcrpiZmNs0RQvk428Zqs3bhszsIkmDkh50928nfq9u54LOqls+zNTcLrugD0u6YsLXl7didXbczBZIUuvjiS7PJxMzm6XxhH/Z3V9rhWt5LhVBblfETM7tsgv6e5IWm9kiM7tQ0i8lbS95DkXbLume1uf3SHqzi3PJxMZvkfKipAPu/sKEb9XuXCqE3K6AmZ7bpS8sMrPbJP2HpAskbXH3fyt1AjmY2SuSVmh857bjkn4j6Q1Jf5Z0pcZ327vb3c9vLlWKmf1C0l5JH0oaa4Wf0PhrjbU6lyoht7tvpuc2K0UBIBE0RQEgERR0AEgEBR0AEkFBB4BEUNABIBEUdABIBAUdABJBQQeARPw/djJzJlOhSG0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a28cdd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "img = gen_img();\n",
    "img2 = gen_noise_add(img, flip=noise_flip)\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(img)\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(img2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_bitstring(img):\n",
    "    v = []\n",
    "    bs = sdmlib.Bitstring.init_ones(1000)\n",
    "    for py in xrange(height):\n",
    "        for px in xrange(width):\n",
    "            pixel = img.getpixel((px, py))\n",
    "            value = sum([int(x/255+0.5) for x in pixel[:3]])//3\n",
    "            assert value == 0 or value == 1\n",
    "            idx = px+width*py\n",
    "            assert idx >= 0 and idx < 1000, 'Ops {} {} {}'.format(x, y, idx)\n",
    "            bs.set_bit(idx, value)\n",
    "            v.append(value)\n",
    "    v2 = [bs.get_bit(i) for i in xrange(height*width)]\n",
    "    assert v == v2\n",
    "    return bs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_img(bs):\n",
    "    img = Image.new('RGBA', (30, 30), (255, 255, 255))\n",
    "    draw = ImageDraw.Draw(img)\n",
    "    for py in xrange(height):\n",
    "        for px in xrange(width):\n",
    "            idx = px+width*py\n",
    "            assert idx >= 0 and idx < 1000, 'Ops {} {} {}'.format(x, y, idx)\n",
    "            x = 255*bs.get_bit(idx)\n",
    "            draw.point((px, py), fill=(x, x, x))\n",
    "    return img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bits = 1000\n",
    "sample = 1000000\n",
    "scanner_type = sdmlib.SDM_SCANNER_OPENCL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "address_space = sdmlib.AddressSpace.init_from_b64_file('sdm-letters.as')\n",
    "counter = sdmlib.Counter.load_file('sdm-letters')\n",
    "sdm = sdmlib.SDM(address_space, counter, 451, scanner_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read(letter, n=6, p=0.25):\n",
    "    n = 20\n",
    "    cols = 7\n",
    "    rows = n//cols + 1\n",
    "    plt.figure(figsize=(20,10))\n",
    "\n",
    "    img = gen_img(letter=letter);\n",
    "    img2 = gen_noise_add(img, p=p, flip=noise_flip)\n",
    "    plt.subplot(rows, cols, 1)\n",
    "    plt.imshow(img2)\n",
    "\n",
    "    for i in xrange(n):\n",
    "        bs2 = to_bitstring(img2)\n",
    "        bs3 = sdm.read(bs2)\n",
    "        if bs3 == bs2:\n",
    "            break\n",
    "        img3 = to_img(bs3)\n",
    "        plt.subplot(rows, cols, i+2)\n",
    "        plt.imshow(img3)\n",
    "        img2 = img3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIEAAACqCAYAAAA6El8nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGQ1JREFUeJzt3XusnWWd6PHfj8sELxDLoRKudlILiOgpyYYA42VOuKRHjUhKYIhOGjXTCY4JJiSKnFhBY0Ki43jiOY5AJAXjAU9SJhKvOEjEwWrcQwgXYaAgEy6lLXoMokEFnvNHF9q1393utddel+fd7+eTNOzn6bu6nr3W17Xjk7WenaWUAAAAAGB522/aCwAAAABg/GwCAQAAAHSATSAAAACADrAJBAAAANABNoEAAAAAOsAmEAAAAEAH2AQCAAAA6ACbQAAAAAAdsKRNoMxcl5n/kZnbMvOyUS0KhqFHaqFFaqFFaqJHaqFFaqFFpiFLKcPdMHP/iHgoIs6OiCci4mcRcVEp5ed7u81hhx1WVq1aNdT9sfw99thj8cwzz+Qwt11sj1pkXybZYoQe2bdhe9Qio+bnNLXwc5qa+DlNLQZt8YAl3MepEbGtlPJoRERm3hQR50bEXqNdtWpVzM7OLuEuWc5mZmaWcvNF9ahF9mWSLUbokX1bQo9aZKT8nKYWfk5TEz+nqcWgLS7l42BHRcTje4yf6M31ycyNmTmbmbO7du1awt3BPi3YoxaZEK+N1EKL1MTPaWrhtZFaaJGpGPvB0KWUa0opM6WUmZUrV4777mCvtEhN9EgttEgttEhN9EgttMioLWUT6MmIOGaP8dG9OZgGPVILLVILLVITPVILLVILLTIVS9kE+llErMnMv8zMv4iIv4mIW0azLFg0PVILLVILLVITPVILLVILLTIVQx8MXUp5ITM/HBHfi4j9I+K6Usr9I1sZLIIeqYUWqYUWqYkeqYUWqYUWmZal/HawKKV8OyK+PaK1wJLokVpokVpokZrokVpokVpokWkY+8HQAAAAAEyfTSAAAACADrAJBAAAANABNoEAAAAAOsAmEAAAAEAH2AQCAAAA6ACbQAAAAAAdYBMIAAAAoANsAgEAAAB0gE0gAAAAgA6wCQQAAADQATaBAAAAADrAJhAAAABAB9gEAgAAAOiAA5Zy48x8LCJ+ExEvRsQLpZSZUSwKhqFHaqFFaqFFaqJHaqFFaqFFpmFJm0A9/62U8swI/h0YBT1SCy1SCy1SEz1SCy1SCy0yUT4OBgAAANABS90EKhFxa2b+e2ZunO+CzNyYmbOZObtr164l3h3s0z571CIT5LWRWmiRmvg5TS28NlILLTJxS/042FtKKU9m5msj4vuZ+WAp5Y49LyilXBMR10REzMzMlCXeX7V+8pOf9I1vuummxjV33nln33jbtm2Na5577rnG3Ctf+cq+8bHHHtu45tRTT+0br1+/vnHNunXr+sb77bfs3gi2zx7b1mJmLnhNKc1vY5Db1Wa+76PlWv/a2MaOBrEMW1tIq1pcLt11sLNBtebntBaXPa+NU6DHeWlxCrre4pJ2AUopT/b+uzMi/iUiTt33LWB89EgttEgttEhN9EgttEgttMg0DL0JlJmvysyDX/46Is6JiPtGtTBYDD1SCy1SCy1SEz1SCy1SCy0yLUv5ONjhEfEvvbeEHRAR/6eU8t2RrAoWT4/UQovUQovURI/UQovUQotMxdCbQKWURyPiv45wLTA0PVILLVILLVITPVILLVILLTItSz0YuhPuv//+vvGHPvShxjV33HFHY25Unn322b7xffc13yU4d+66665rXHPCCSf0jb/4xS82rjnrrLOGWSIjMMxBa8vlcLb5vo+uH9g2Sculo0EMe+A6LMao/jelRZbKax410SO16HqLy+7XQwEAAADQZBMIAAAAoANsAgEAAAB0gE0gAAAAgA5wMPQcW7Zsacy9973v7Rv//ve/n9RyRurBBx/sG5999tmNa6688srG3KZNm8a2pq7q0kG80DYOKge6pOsHpFIXPVKL5dyidwIBAAAAdIBNIAAAAIAOsAkEAAAA0AGdPxPotttu6xtfeOGFjWtefPHFBf+dAw88sG+8YcOGxjXr16/vG5988smNa1asWNGYe/rpp/vGDz/8cOOar371q33jr3/9641rnn/++cbcXJ/85CcbcwcddFDf+KMf/eiC/w4slXNZ2me+52fu8zjINeM07ftnMkb1WjHONryedYMWqYkeqYUWp8s7gQAAAAA6wCYQAAAAQAfYBAIAAADogAU3gTLzuszcmZn37TF3aGZ+PzMf7v23eZANjIEeqYUWqYUWqYkeqYUWqYUWqc0gB0Nvjoj/FRE37DF3WUTcVkq5KjMv640/Nvrljdazzz7bmHvf+97XNx7kEOiTTjqpMbdly5a+8XHHHbfI1e3dscceu89xRMSZZ57ZN/7EJz7RuOa8887rG997770D3f/ll1/eNz7nnHMa16xdu3agf2sENscy6XEU2nBg2TI+dHdzaHFBgzRaW8e1rWcAm0OLY9HCFmqwOfQ4clocyubQ4ljocdE2hxbHQovDWfCdQKWUOyLiV3Omz42I63tfXx8R7xnxumBeeqQWWqQWWqQmeqQWWqQWWqQ2w54JdHgpZXvv66cj4vC9XZiZGzNzNjNnd+3aNeTdwT4N1KMWmQCvjdRCi9TEz2lq4bWRWmiRqVnywdBl93uw9vo+rFLKNaWUmVLKzMqVK5d6d7BP++pRi0yS10ZqoUVq4uc0tfDaSC20yKQNcibQfHZk5hGllO2ZeURE7Bzlosblc5/7XGPu6aefXvB2a9as6Rv/6Ec/alzzmte8ZviFjcHq1asbc1u3bu0bz8zMNK558MEHG3Nzz0n62MeaH1f93ve+t9gljlIre+yK+T6rO+w5QXNvV+HngLXYQhV2NApapCZ6pBZapBZaZGqGfSfQLRGxoff1hoj4xmiWA0PRI7XQIrXQIjXRI7XQIrXQIlMzyK+IvzEitkbE8Zn5RGZ+MCKuioizM/PhiDirN4ax0yO10CK10CI10SO10CK10CK1WfDjYKWUi/byV2fuZR7GRo/UQovUQovURI/UQovUQovUZskHQwMAAABQv2EPhq7eCy+80Ji7+uqrh/q3Nm/e3Deu7RDoQb3qVa/qG99www2Na0477bTG3EsvvdQ3vvXWWxvXPPLII33j+Q6mhpfNPYh32IOiAQAAGJx3AgEAAAB0gE0gAAAAgA6wCQQAAADQAcv2TKCtW7c25nbu3Lng7U4//fTG3BlnnDGSNdXmlFNOacy95S1vaczdcccdC/5b3/zmN/vGl1xyyfALo3PmnhEEAADA6HknEAAAAEAH2AQCAAAA6ACbQAAAAAAdYBMIAAAAoAOW7cHQd95551C3O+ecc0a8kna54IILGnODHAw9Ozs7juUwhwOUAQAAGJZ3AgEAAAB0gE0gAAAAgA6wCQQAAADQAQueCZSZ10XEuyJiZynlpN7cFRHxdxGxq3fZ5aWUb49rkcN4+OGHh7rdm9/85hGvpF3e8IY3DHW7X/ziFyNeyfza2iPLjxbHJzOHul1Xz8zS4nCG7WxYXelTj4unxfHQ4nD0OHpaHI4Wx2eQdwJtjoh188z/Uyllbe+PYJmUzaFH6rA5tEgdNocWqcfm0CN12BxapA6bQ4tUZMFNoFLKHRHxqwmsBRakR2qhRWqhRWqiR2qhRWqhRWqzlDOBPpyZ92TmdZm5Ym8XZebGzJzNzNldu3bt7TJYqgV71CIT4rWRWmiRmvg5TS28NlILLTIVw24C/XNErI6ItRGxPSL+cW8XllKuKaXMlFJmVq5cOeTdwT4N1KMWmQCvjdRCi9TEz2lq4bWRWmiRqVnwYOj5lFJ2vPx1Zl4bEd8c2YpG5Ne//vVQtzv00ENHvJJ2OfLII4e63W9/+9sRr2RwbeiRbtDiwiZ9yF9XabFJe9Ojx35anB4tNulxOrTYpMXJGeqdQJl5xB7D8yLivtEsBxZPj9RCi9RCi9REj9RCi9RCi0zTIL8i/saI+OuIOCwzn4iIT0bEX2fm2ogoEfFYRPz9GNcIf6JHaqFFaqFFaqJHaqFFaqFFarPgJlAp5aJ5pr8yhrXAgvRILbRILbRITfRILbRILbRIbZby28EAAAAAaImhDoZug1LKtJfQSsMeyHXIIYeMeCXUaG4f/nfGYs3XjIMAmYS57U26O6+f1EKL1ESP1KJLLXonEAAAAEAH2AQCAAAA6ACbQAAAAAAdsGzPBFqxYsVQt/vlL3854pW0y7Df/8qVK0e8Emq0nD8by/ToimkYZXddOkeApRtnH1pksfRILbQ4Od4JBAAAANABNoEAAAAAOsAmEAAAAEAH2AQCAAAA6IBlezD0mjVrhrrd3Xff3Zhbv379UpfTGvN9/4N44xvfOOKVALXr+qF61GPaLU77/qnHtFuY9v1Tl2n3MO37px7TbmHa918b7wQCAAAA6ACbQAAAAAAdYBMIAAAAoAMWPBMoM4+JiBsi4vCIKBFxTSnlf2bmoRHx9YhYFRGPRcQFpZT/N76lLs4ZZ5wx1O2++93vNuY+/elPL3U5rTHf9z+IdevWjXglTW1tcZQys2+8nD/fOvd7nc80v389dqvHmmlRi7XQohZrokc91kKLWqzNIO8EeiEiLi2lnBgRp0XEP2TmiRFxWUTcVkpZExG39cYwTlqkJnqkFlqkFlqkJnqkFlqkKgtuApVStpdS7up9/ZuIeCAijoqIcyPi+t5l10fEe8a1SIjQInXRI7XQIrXQIjXRI7XQIrVZ1JlAmbkqIk6OiJ9GxOGllO29v3o6dr+9bb7bbMzM2cyc3bVr1xKWCn+mRWqiR2qhRWqhRWqiR2qhRWow8CZQZr46IrZExEdKKc/u+Xdl94f65v1gXynlmlLKTCllZuXKlUtaLERokbrokVpokVpokZrokVpokVoseDB0RERmHhi7g/1aKeXm3vSOzDyilLI9M4+IiJ3jWuQw5jsY+rWvfW1jbufO/mXPzs42rvnhD3/YN37729++xNUtzmc+85m+8amnntq45swzz+wb77ffwvt799xzT2PuW9/61oK3W716dWNuvjWNQxtbZGGDHAJdIz1SCy1SCy1SEz1SCy1SkwV3CnL3/zv7SkQ8UEr5/B5/dUtEbOh9vSEivjH65cGfaZGa6JFaaJFaaJGa6JFaaJHaDPJOoL+KiL+NiHsz8+7e3OURcVVE/N/M/GBE/GdEXDCeJcKfaJGa6JFaaJFaaJGa6JFaaJGqLLgJVEr5t4jY22c1ztzLPIycFqmJHqmFFqmFFqmJHqmFFqnNQGcCtdEBBzS/tYsvvrgxd+WVVy74b33gAx/oG893btCKFSsWsbq927FjR2Nu06ZNfeOXXnqpcc3rX//6vvHVV1/duOaUU07pG8/9vvb2b8/1qU99qjG3//77L3g7eNmwZwDtPjMPAACAYSzqV8QDAAAA0E42gQAAAAA6wCYQAAAAQAfYBAIAAADogGV7MPR8Lr300sbctdde2zd+6qmnGtc8+uijfeO3vvWtjWtuvvnmvvFxxx03zBLj+eefb8y9//3v7xtv2bKlcc22bdv6xmeddVbjmte97nV948cee2ygNb373e/uG1900UUD3Y5uGvbQZwAAAMbLO4EAAAAAOsAmEAAAAEAH2AQCAAAA6ACbQAAAAAAd0KmDoQ8++ODG3I033tg3nu9A5T/+8Y994/vvv79xzZve9Ka+8YYNGxrXnH/++X3jtWvXNq458sgjG3NXXXVV3/id73xn45qNGzf2jZ955pnGNYMcBH388cc35uY+Rg7+rYfngpq0ocdSyrSXQMUm2bAW2Zdxtqg9FkuP1EKLo+GdQAAAAAAdYBMIAAAAoAMW3ATKzGMy8/bM/Hlm3p+Zl/Tmr8jMJzPz7t6fd4x/uXSZFqmJHqmFFqmFFqmJHqmFFqnNIGcCvRARl5ZS7srMgyPi3zPz+72/+6dSyufGt7zxe9vb3tY3vvnmmxvXXHjhhX3j3/3ud41r/vCHP/SNr7322sY1883V5qGHHmrMffazn+0bb9q0qXHNhM5RWDYtzveZ0zacp0KfZdNjlwz7v7PKPyeuxTm8nk6NFueoscW5a6r89W0p9DiHHqdGi3NocboW3AQqpWyPiO29r3+TmQ9ExFHjXhjMpUVqokdqoUVqoUVqokdqoUVqs6gzgTJzVUScHBE/7U19ODPvyczrMnPFXm6zMTNnM3N2165dS1osvEyL1ESP1EKL1EKL1ESP1EKL1GDgTaDMfHVEbImIj5RSno2If46I1RGxNnbvbP7jfLcrpVxTSpkppcysXLlyBEum67RITfRILbRILbRITfRILbRILQbaBMrMA2N3sF8rpdwcEVFK2VFKebGU8lJEXBsRp45vmbCbFqmJHqmFFqmFFqmJHqmFFqnJgmcC5e4Tkr4SEQ+UUj6/x/wRvc83RkScFxH3jWeJk/Wud72rMXfXXXf1jS+++OLGNbfffvvY1jQqBx10UN/4+eefb1wz3wFYV1xxRd/4xz/+ceOaL33pS33j1atXD7HCfVvuLc597Gs8MG1c2njw2nLq0UHl/drWoxbboW1dDUOL7dCFFiP02BZd6FGL7dCFFl82yG8H+6uI+NuIuDcz7+7NXR4RF2Xm2ogoEfFYRPz9WFYIf6ZFaqJHaqFFaqFFaqJHaqFFqjLIbwf7t4iYb3vv26NfDuydFqmJHqmFFqmFFqmJHqmFFqnNon47GAAAAADtNMjHwTrv+OOP7xv/4Ac/aFyzdevWvvGNN97YuObOO+/sGz/yyCONa5577rnG3Cte8Yq+8dFHH9245uSTT+4br1u3rnHN+vXr+8bf+c53Gtdccskljbmnnnqqb3zrrbc2rjnppJP6xo8//njjmsMOO6wxx94N8pnbtn4ut0ufuW2rNp5RpavlqQ0taq8btEhN9EgttNg+3gkEAAAA0AE2gQAAAAA6wCYQAAAAQAfYBAIAAADoAAdDj8jpp5++z3GNzj///MbcOeec05j7whe+0Df+8pe/3LjmhBNO6Bs7BHo8BjnUzMFnjIOuqIUWqYUWqYkeqYUW6+edQAAAAAAdYBMIAAAAoANsAgEAAAB0gDOB6HPIIYc05jZt2tQ3/vjHP964ZseOHWNbEwAAALB03gkEAAAA0AE2gQAAAAA6wCYQAAAAQAfYBAIAAADogCylTO7OMndFxH9GxGER8czE7nh02rjuNq35daWUlZO4oz1ajGjXY/Qyax6vibUY0frXxjauOaJd657Ga2ObHp89tXHdbVqzFgfXxnW3ac1+Tg+ujWuOaNe6vTYOpo1rjmjXugdqcaKbQH+608zZUsrMxO94idq47jauedLa+BhZ8/LUxseojWuOaO+6J6Wtj08b193GNU9SWx+fNq67jWuetDY+Rm1cc0R71z0pbXx82rjmiPaue198HAwAAACgA2wCAQAAAHTAtDaBrpnS/S5VG9fdxjVPWhsfI2tentr4GLVxzRHtXfektPXxaeO627jmSWrr49PGdbdxzZPWxseojWuOaO+6J6WNj08b1xzR3nXv1VTOBAIAAABgsnwcDAAAAKADbAIBAAAAdMDEN4Eyc11m/kdmbsvMyyZ9/4PIzOsyc2dm3rfH3KGZ+f3MfLj33xXTXONcmXlMZt6emT/PzPsz85LefNXrnqY2tBihxy7Q4vhocfHa0KMWu6ENLUbosQu0OD5aXLw29KjFuk10Eygz94+I/x0R/z0iToyIizLzxEmuYUCbI2LdnLnLIuK2UsqaiLitN67JCxFxaSnlxIg4LSL+offY1r7uqWhRixF6XNa0OHZaXIQW9bg5tListajFCD0ua1ocOy0uQot63BxarNak3wl0akRsK6U8Wkr5Q0TcFBHnTngNCyql3BERv5ozfW5EXN/7+vqIeM9EF7WAUsr2Uspdva9/ExEPRMRRUfm6p6gVLUbosQO0OEZaXLRW9KjFTmhFixF67AAtjpEWF60VPWqxbpPeBDoqIh7fY/xEb64NDi+lbO99/XREHD7NxexLZq6KiJMj4qfRonVPWJtbjGjR86rHBWlxQrQ4kDb32JrnVIsDaXOLES16XvW4IC1OiBYH0uYeW/OcLvcWHQw9hFJKiYgy7XXMJzNfHRFbIuIjpZRn9/y7mtfN8Gp+XvXYLTU/p1rslpqfUy12T83Pqx67pebnVIvdUvNz2oUWJ70J9GREHLPH+OjeXBvsyMwjIiJ6/9055fU0ZOaBsTvYr5VSbu5NV7/uKWlzixEteF71ODAtjpkWF6XNPVb/nGpxUdrcYkQLnlc9DkyLY6bFRWlzj9U/p11pcdKbQD+LiDWZ+ZeZ+RcR8TcRccuE1zCsWyJiQ+/rDRHxjSmupSEzMyK+EhEPlFI+v8dfVb3uKWpzixGVP696XBQtjpEWF63NPVb9nGpx0drcYkTlz6seF0WLY6TFRWtzj1U/p51qsZQy0T8R8Y6IeCgiHomI/zHp+x9wjTdGxPaI+GPs/pzlByPiv8Tu08Afjoh/jYhDp73OOWt+S+x+a9o9EXF37887al/3lB+z6lvsrVOPy/yPFse6Zi0u/jGrvkctduNPG1rsrVOPy/yPFse6Zi0u/jGrvkct1v0ne98wAAAAAMuYg6EBAAAAOsAmEAAAAEAH2AQCAAAA6ACbQAAAAAAdYBMIAAAAoANsAgEAAAB0gE0gAAAAgA74/0GvEu6ZswoSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a34a590>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "read('Q', p=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
